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Abstract 

In many mathematical models of physical phenomenons and engineering fields, such as electrical 
circuits or mechanical multibody systems, which generate the differential algebraic equations 
(DAEs) systems naturally. In general, the feature of DAEs is a sparse large scale system of fully 
nonlinear and high index. To make use of its sparsity, this paper provides a simple and efficient 
algorithm for computing the large scale DAEs system. We exploit the shortest augmenting path 
algorithm for finding maximum value transversal (MVT) as well as block triangular forms (BTF). 
We also present the extended signature matrix method with the block fixed point iteration and 
its complexity results. Furthermore, a range of nontrivial problems are demonstrated by our 
algorithm. 

Key words: Differential algebraic equations, sparsity, shortest augmenting path, block 
triangular forms, structural analysis 


1. Introduction 


The problem of differential algebraic equations (DAEs) system solving is fundamental in 
modelling many equation-based models of physical phenomenons and engineering fields, such 
as electric circuits [28i, 29|, mechanical systems E^l . spacecraft dynamics llzil , chemical engi¬ 


neering J32], and many other areas. Generally, DAEs can be produced very large scale system 
of fully nonlinear and higher index in practice. However, most of the algorithms treat the low 
index case or consider solutions of linear systems 114, JJ, 19, 3(3]. The index is a notion used 


in the theory of DAEs for measuring the distance from a DAE to its related ordinary differential 
equation (ODE). It is well known that it is direct numerical computations difficult to solve a high 
index DAE. In particular, it may only solve some special classes of DAEs by the direct numerical 
solution floH lal. 

Index reduction techniques can be used to remedy the drawback of direct numerical compu¬ 
tation |3|]. Pantelides’ method 0 gives a systematic way to reduce high index systems of DAEs 
to lower index one, by selectively adding differentiated forms of the equations already present in 
the system. In j^l. Ding et al. developed the weighted bipartite algorithm based on the minimally 
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structurally singular subset, which is similar to the Pantelides’ method. However, the algorithms 
can succeed yet not correctly in some instances 1271. Campbell’s derivative array 0] needs to 
be computationally expensive especially for computing the singular value decomposition of the 
Jacobian of the derivative array equations using nonlinear singular least squares methods. Sig¬ 
nature matrix method (also called 2-method) 12111 is based on solving an assignment problem, 
which can be formulated as an integer linear programming problem. In [211], Pryce proved that 2- 
method is equivalent to the famous method of Pantelides’ algorithm, and in particular computes 
the same structural index. However, the nice feature of 2 method is a simple and straightforward 
method for analyzing the structure of DAEs of any order, not just first order. 

In particular, large scale and high index DAEs with fully nonlinear systems are now routine 
that such models are built using interactive design systems based on the Modelica language 
[7, Il2ll . In addition, the sparsity pattern of DAEs can arise in most actual applications 112 lL 
l24ll . In ifTltl . Frenkel et al. gave a survey on appropriate matching algorithms based on the 
augmenting paths and push-relabel algorithm by translating Modelica models for large scale 
systems of DAEs. More recently, Wu et al. [[33] generalized the 2-method to the square and 
/-dominated partial differential equations (PDEs) systems. Pryce et al. [23] generalized the 2- 
method for constructing a block triangular form (BTF) of the DAEs and exploiting to solve it 
efficiently in a block-wise manner. In & Tang et al. proposed the block fixed-point iteration 
with parameter method for DAEs based on its block upper triangular structure. However, the 
essential task is to solve the linear assignment problem for finding a maximum value transversal 
(MVT), which is a large part of the cost for index reduction of DAEs solving. Pryce et al. 
mentioned only in their work using Cao’s Matlab implementation 0] of the shortest augmenting 
path algorithm of Jonker and Volgenant in 0. We focus on solving in the large scale and high 
index cases in order to provide the shortest augmenting path algorithm for finding an MVT and 
an extended signature matrix method. The problem is also closely related to computing the block 
triangular form of a sparse matrix and linear assignment problems over integer. 

Our approach is based on signature matrix method and modified Dijkstra’s shortest path 
method. Our fundamental tool is the block triangularization of a sparse matrix; we exploit re¬ 
cent advances in linear assignment problem solving, which is equivalent to finding a maximum 
weight perfect matching in a bipartite graph of signature matrix in 2-method, and we adapt the 
block fixed-point iteration with parameter for the canonical offsets techniques. Currently, we 
are working on the theoretical foundation and implementation of these methods on Maple plat¬ 
form. Another direction for future work is to exploit the fact that our algorithms are expressed in 
OpenModelica solvers. 

The rest of this paper is organized as follows. The next section introduces our purpose and the 
shortest augmenting paths based algorithms, and presents an improved algorithm for the block 
triangularization for DAEs system. Section 3 describes the extended signature matrix method for 
the structural analysis of large scale DAEs system and gives its complexity results. The following 
section shows our algorithm for an actual industrial example and some experimental results. The 
final section makes conclusions. 
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2. Preliminary results 

2.1. Purpose 

We consider an input DAEs system in n dependent variables Xj = Xj(t ) with 1 a scalar inde¬ 
pendent variable, of the general form 


fi(t, the Xj and derivatives of them) = 0, 1 < i, j < n. (1) 

The f are assumed suitably smooth, and the derivatives of xj are arbitrary order. In general, 
signature matrix method is an effective preprocessing algorithm for the small and middle scale 
DAEs system. First, it needs to form the n x n signature matrix I = of the DAEs, where 

( highest order of derivative to which the variable Xj appears in equation f, 
or — oo if the variable does not occur. 

Then, taken the analysis procedure of E as a linear assignment problem is to seek the offsets of 
the DAEs, that is, the number of differentiations of f. It can be formulated by the following 
primal problem: 


Maximize z — E (A/ft/, 

(>J)eS 

sub ject to E &j = 1 V i = 1,2, ■ • • , n, 

MhfleS) 

E fij — 1 V 7 = 1,2, ■■■,*, 

{0,1} V O', j) e S. 

Note that the state variable tf } only be defined over the sparsity pattern of the problem 


( 2 ) 


S = sparse (I.) = !(i, j)|<x, ; > -oo). (3) 

It can be also defined on an undirected bipartite graph, in which case an assignment is a perfect 
matching. Given a bipartite graph G(E) — (F,X,e), where F is the set vertices corresponding to 
the rows of E, X is the set vertices corresponding to the columns of E, and e is the set of edges 
corresponding to the non-negatively infinite in E, f’| = |X| = n, | • | denotes the cardinality of a 
set. In this paper, our goal is to handle large scale systems with n involving thousands and even 
more. 


2.2. Block triangularization for DAEs system 

As we have encountered with the increasingly large problems, an important preprocessing 
manipulation is the block triangularization of the system, which allows to split the overall system 
into subsystems which can be solved in the sequence or parallelization jg], Pothen et al. j20f| 
described implementations of algorithms to compute the block triangular form by permuting the 
rows and columns of a rectangular or square, unsymmetric sparse matrix. It is equivalent to 
computing a canonical decomposition of bipartite graphs known as the Dulmage-Mendelsohn 
decomposition. Considering the index reduction for DAEs system, the block triangularization 
algorithm can be directly performed on the signature matrix E. The block triangular form of E 
can be generated by the block triangularization of the incidence matrix of the sparse pattern S 
for a given DAEs system. The incidence matrix is defined as follows: 
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Definition 2.1. Let A = [ a,j] be an incidence matrix of the sparse pattern S for a given DAEs 
system, whose rows and columns represent the F and X as above, and element a,j is 1 if(i, j) e S 
and is 0 otherwise. 

The associated bipartite graph G{A) for the incidence matrix A is defined as follows: 

Definition 2.2. Let G(A) — ( F,X,E ) be the associated undirected bipartite graph for the in¬ 
cidence matrix A of a given DAEs system, where F and X are defined as above, one for each 
equation f, and other for each variable Xj respectively; and ( f,Xj) belongs to E if and only if 
a jj-t~ 0. 

Definition 2.3. A subset M of E is called a matching if a vertex in G(A ) is incident to at most 
one edge in M. A matching M is called maximal, if no other matching M' D M exists, which 
is called maximum if\M\ > \M'\ for every matching M'. Furthermore, if\M\ — |F| = \X\, M is 
called a perfect matching. 

Lemma 2.1. iiTH/i A bipartite graph G(A) with vertex sets F and X contains a perfect matching 
from F to X if and only if it satisfies Hall’s condition 

|r(/)| > 1/1 for every f C F, 

where T(f) — {xj e X, (f, xj) 6 E for some f e F) c X. 

A walk is a sequence of nodes vo, Vi,--- , v„_i e F!JX suc h that (v,-,v,+i) e E for i = 
0,1, ■ • ■ , n - 2. Noted that edges or nodes can be repeated in a walk. An alternating walk is a 
walk with alternate edges in a matching M. An alternating path is an alternating walk with no 
repeated nodes. With respect to M, we can define the following sets: 

VF — {F-nodes reachable by alternating path from some unmatched F-node) 

HF = IF-nodes reachable by alternating path from some unmatched X-node) 

SF = F\(VF\JHF) 

VX = {X-nodes reachable by alternating path from some unmatched F-node) 

HX = IX-nodes reachable by alternating path from some unmatched X-node! 

SX = X\(VX{JHX) 

If there exists a perfect matching for G(A), the directed graph G,/(A) of G(A ) is made by 
the rule: let non-matching edges in G(A) be directed from X-nodes to F-nodes, matching edges 
shrunk into single nodes, and the nodes identified with F. In order to demonstrate the bene¬ 
fits of parallel processing for index reduction, we present an improved algorithm for the block 
triangularization of the incidence matrix A. The main steps are as follows: 

step 1 Find the diagonal sub-blocks of A from the connected components of G(A). Assume it 
generates p diagonal submatrices A\, A 2 ,..., A p . It is helpful to use parallel computing for 
each A,-. 

A 1 

A 2 0 

0 A p ~ 1 

Ap] 
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step 2 If A has a perfect matching, then it only has the square A s , otherwise the underdetermined 
A/, and overdetermined A,, will be present. For each A,, find a maximum matching M, in 
the associated graph G(A,). 

step 3 With respect to M,, partition F, into the sets VF,. SFj, IIF, using the rule of Dulmage- 
Mendelsohn decomposition; partition X, into the sets VX,, SXj, IIX, similarly, where F, 
and X, are the block form of F and X respectively, 'x' denotes a possibly nonzero matrix 
of appropriate dimensions. 

HXi SXi VX, 

A;, X X 

0 A h x 

0 0 A,, . 


Step. 3 HF, 

Ai —* S Fj 

VFj 


step 4 Using Tarjan’s depth-first search algorithm, find the block upper triangular form of the 
square submatrix A, s (A,; and A,- (i are not square) by finding strong components in the 
associated directed graph GfiAj,). Assume it produces A it j, A, s g, ..., A, sj „ ; . 

sx Uni 

x 
x 

Ai„ ni - 



SX,j 

SXi, 2 

S Fj\ 

A;„i 

X 

StepA S Fjg 

A, s —* 

0 

Ai s , 2 

SF un , 

0 

0 


Remark 2.2. Some problems may only require one step of the computation in step 1, one at 
a time; others can take advantage of the further decomposition. In addition, if Mi is perfect 
matching(Aj contains at least one transversal), then HX„VXj,HF) and VFj are empty, that is, 
the step 3 can be skipped and A, = Aj s . Furthermore, our algorithm can be generalized to the 
underdetermined and overdetermined signature matrix of sparse pattern S case. 

Remark 2.3. (a) The following are equivalent: 

(i) A is structurally nonsingular (contains at least one transversal). 

(ii) A has the Hall’s condition from Lemma \2J\ 

(b) The following are equivalent: 

(i) A is irreducible. 

(ii) A has the strong Hall’s condition, which satisfies for all proper subsets |r(/)| > j/| + 1. 

(iii) Every element of A is on a transversal, and the bipartite graph of A is connected. 

Theorem 2.4. The above algorithm works correctly as specified and its complexity as follows: 

(a) If A is irreducible, then it only performances in step 1, that is, 0(n + t ) operations, 

(b) If A is reducible, then it needs the whole steps, that is, 0(pN <X>) operations, 
where N and <I> are the constant, as defined below. 
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Proof. Correctness of the algorithm follows from the basic idea. The essential tasks of above 
algorithm are to permute the given matrix to block triangular form. A permuted block form of the 
DAEs system is generated from forming permutations F and X of the equations and variables, 
which denote the block form (Fi, F 2 , • • • , F p ) and (X \, X 2 , • ■ • , X p ) with n\ + 112 + • • ■ + n p = n, 
where |F,j = |A,j = i - 1,2, • • ■ , p. 

From the description of algorithm, we observe that there are four major steps on time com¬ 
plexity. In step 1, we can use a version of recursive depth-first search to find the connected 
components in the run time 0(n + r), where r is the number of nonzero entries in A (r = 
ti + T 2 + • ■ • + T p ). Suppose n, = N, t, — Q, for each i, that it, n = p ■ N, r = p • <f>, a 
0{N<S>) algorithm for each A,-(i = 1,2, ■ • ■ ,p) is considered in step 2, which can be parallel com¬ 
puting. It terminates with a maximum matching in the graph. In step 3, it only needs to check 
the Mi form, which costs in little time. In step 4, Tarjan’s depth-first search algorithm to find 
the strong components of a directed graph is in linear time of its edges, that is, the bound to the 
optimal 0(N + <I>). In all, the worst-case complexity of the algorithm becomes 0(pN<\>). □ 

The E of sparse pattern S for a given DAEs system can be block triangulated by the trian- 
gularization algorithm mentioned as above. Therefore, the E with block triangular form is as 
follows: 


E 


1.1 


E = 


St,2 
^ 2,2 



(4) 


E 


p,p 


where the elements in the blanks of E are - 00 . In general, there are three possible cases for E of 
sparse pattern S as follows, see Figures 1, 2 and 3. 



Figure 1: irreducible E Figure 2: E with only diagonal Figure 3: E with block triangu- 

blocks lar form 

From Figures 1, 2 and 3, we have the following observations. 

• For the E of sparse pattern S , there may exist some strong coupling terms with the block 
triangular form in Figure 1 other than Figures 2, 3. It cannot be permuted to a nontrivial 
BTF (one with p > 1), whose matrix is said to be irreducible, although there exists the 
maximum number of nonzeros on the diagonal of E. If the E can be permuted to the BTF 


1 /V. <I> arc defined by the same way for the rest of this paper. 
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Q with p > 1 in Figures 2, 3 is called to be reducible. Therefore, if the X contains one 
transversal, that is, S for the X contains one transversal, then the X, y is irreducible for 
i, j = 1,2, ■ • • , p. If the X for the given DAEs system contains some transversal, the BTF 
of X can be found by the triangularization algorithm as above, and its diagonal blocks are 
square and irreducible. 

• For Figure 2, each element of each diagonal block is on a transversal from Remark[273jb)(iii) 
in the absence of coupling. It can lead to the parallelization for computing the transver¬ 
sal on each block naturally. The X contains a remarkable property, which of difference 
between global offsets of X and local offsets on each block must be the constant. 

• For Figure 3, the X can be permuted to the block upper triangular form. It will be influenced 
by the top right blocks for computing the transversal, since the lower left blocks are blank. 
Moreover, it needs the block fixed point iteration method with parameter for computing 
the offsets of X from top to bottom in sequence. In particular, for all column j if <x l; - in the 
top right blocks are less than the elements of the same columns in the diagonal blocks all 
the time, then it is similar to Figure 2. 

2.3. Shortest augmenting path algorithm 

In this subsection, we present the shortest augmenting path algorithm to find a maximum 
value transversal in the sparse pattern S of signature matrix X = from Q. A transversal of 
S is an n element subset T of S with just one element in each row and column. The DAEs system 
is structurally well-posed if there exists a T, all of whose <t (/ are finite, else structurally ill-posed. 
The problem is to find a maximum value transversal, which makes ||7j| = fti.per (T ij as large as 
possible. Meanwhile, it is equivalent to finding a maximum cardinality matching in a bipartite 
graph whose incidence matrix is the signature matrix, and can be solved by the AlgorithmU] 

Remark 2.5. The complexity of the classical algorithm for finding maximum value transversal 
can be done at Off) operations based on the Hungarian method. For instance, Balinski 10/ pre¬ 
sented a signature methods for the assignment problem, which considers feasible dual solutions 
corresponding to trees in the bipartite graph of row and column nodes. The shortest augmenting 
path algorithm can be reduced to 0(n~ log n) operations by using priority queues. By fully ex¬ 
ploiting monotonically nondecreasing over all augmentations for the minimum distance, it can 
be reduced the average time per augmentation to Ofi). For the sparse case, the whole complex¬ 
ity is more like Off) operations. Moreover, we give a necessary preprocessing to the negative 
of the original X, since the shortest augmenting path algorithm is to compute the minimum-cost 
network flow problem, A parallel version of Algorithmic an speed up the MVTfinding procedure 

a. 

Example 2.1. We present an index-five model of the planar motion of a crane to illustrate the 
Algorithm\ 7] The model is discussed in j^], see Figure 0 It is aimed to determine the horizontal 
velocity u\ of a winch of mass Mi, and the angular velocity u 2 of the winch so that the attached 
load M 2 moves along a prescribed path pft) and p 2 f)- 

The resulting DAEs are given by the Newton’s second law of motion as follows, 

0 = /1 = M 2 x + Tsind, 0 = f 2 — Mfz + tcosO - mg, 

0 = fo = M\d + C\d - u\ - rsind, 0 = fy — Jr + C 2 r + C 3 U 2 — C\t, 

0 = fs - rsind + d - x, 0 = = rcosO - z, 

0 = fi = x- pft), 0 = / s = z — p 2 f). 
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Algorithm 1 (Shortest augmenting paths based algorithms) 

Input: signature matrix X = (cr, ; ). 

Output: either output a maximum value transversal T for X or give a nonexistent error. 
Step 1: Initiation. 

1 : Column reduction. 

2: for j — n to 1 do 

3: find minimum value over rows 

4: if minimum value in a row appears at an unassigned column then 

5: initialize the assignment column 

6 : end if 

7: end for 

8 : Reduction transfer. 

9: for each assigned row i do 

10: transfer from unassigned to assigned rows to a column k 

it: end for 

12 : Augmenting row reduction. 

repeat 

finding augmenting paths starting in an unassigned row i 
if j is unassigned then 

augmentation solution is from the alternating path 


else 

reassigning column j to row i, reduction is transferred 

end if 

until no reduction transfer or augmentation. 

Step 2: Augmentation. 

l: repeat 

2 : Augment solution for each free rows. 

3: construct the auxiliary network and determine from an unassigned row i to an unassigned 

column j 

4: find an alternating path by the modified Dijkstra’s shortest path method, which is used to 

augment the solution 
5: Adjust the signature matrix. 

6 : update column values, and reset row and column assignments along the alternating path 

7: until no unassigned rows. 

Step 3: return the sequence of rows, and the corresponding column indices. 


where (;r, z) is the location of the load, t is the tension in the cable, 6 is the angle of the cable 
with the vertical, J is the moment of inertia of the winch, d is the trolley location and r is the 
cable length, C\ > 0, C? > 0 C 3 > 0, g > 0 and m > 0 are constants. 
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d 



Figure 4: Control of a Crane 


From the definition o/Z in Section |2~71 labeled by equations and dependent variables, is 


fi 

fi 

fi 

z = fi 
/s 
/« 
ft 

fi 


x z d r 

2 

2 

2 

2 

0 0 0 

0 0 

0 

0 


6 T U\ U2 

0 0 
0 0 
0 0 0 

0 0 

0 

0 


where a blank denotes — oo. 

From Section 12.21 we have the block triangularization for © as follows: 


( 6 ) 



U 2 u\ d r 

r 6 

z x 

u 

r0 2 

0 


h 

0 2 

0 0 


fs 

0 0 

0 

0 

k 

0 

0 

0 

h 


0 0 

2 

h 


0 0 

2 

h 



0 

fi 



oJ 


which is similar to the form of Figure\3\ 

From Section IZi] we can get the MVT by Algorithm^: (f\, t), (fi, 6), (f$, u\), f/ 4 , U 2 ), (f$, 
d), (h, r), (fi, x), r/g, z). 
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3. An extended E-method 


We generalize the signature matrix method to the large scale DAEs system with block fixed 
point iteration. The basic idea mainly includes the preprocess of the original system with the 
block triangularization, by adopting the shortest augmenting path algorithm for finding maxi¬ 
mum value transversal and solution of dual problem with block fixed point iteration. Without 
loss of generality, we assume that the DAEs system be structurally well-posed, that is, there 
exists a transversal. The major steps are as follows. 


step 1 Form the n X n signature matrix E = (cr, ; ) of the DAEs system with sparse pattern S. 

step 2 Call the block triangularization algorithm in Section 2.2 for dealing with E, and obtain 
the block form E, 7 for i,j = 1,2, • • • , p. 

step 3 Solve an assignment problem to find the MVT T k (k = 1,2, ■ • • , p) from each block 'Lf-k 
by the shortest augmenting path algorithm separately. 


step 4 Determine the local canonical offsets of the problem, which are the vectors (c^) lsk£p — 
(ci)\<i<„ t , (Ak)\<k< P = (dj)\<j<n k such that dj-Ct > cr ih for all 1 < i, j < n k , and d s -c, = cr, 7 
when (i, j ) e T k . The global offsets (c) = (c/)i<;<„, (d) = (dj)\<j<„ such that dj - Ci > (Tjj, 
for all 1 < i, j < n , and d s - c,- = cr, 7 when (i, j ) e T. It is well known that the difference 
between global and local offsets is the constant on each block when the irreducible fine 
BTF is of the Jacobian sparsity pattern [23]. 


This problem can be formulated as the dual of (0) in the variables c = (ci, C 2 , • • • , c n ) and 
d = (d\,d 2 , • ■ • , d n ), the dual is: 


Minimize 
sub ject to 


Z-Yjdj -YjCi, 

j 

dj - Ci > crjj for all ( i , j ) e S, 
c, > 0 for all i. 


(7) 


In order to compute the local canonical offsets of the problem, we apply the fixed point 
iteration method with parameter to process each irreducible diagonal matrix in block upper 
triangulated E from top to bottom in sequence. The theory of block fixed point iteration 
with parameter method is described in detail by the companion paper M3 ill . 


step 5 To verify the success of the index reduction, we need to check whether the n x n system 
Jacobian matrix J is nonsingular, where 


df 


j _ I d((dj—Ci)th derivative of Xj) 


o 


if this derivative is present in f, 
otherwise. 


In this paper the structural index is then defined as: 


0 for all di > 0 

v = max,- Ci + < 

11 for some dj — 0. 

step 6 Choose a consistent point. If J is nonsingular at that point, then the solution can be com¬ 
puted with Taylor series or numerical homotopy continuation techniques in a neighborhood 
of that point. 
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It is well known that if Tj. is the MVT of X« with BTF (J4]», then T = IJ 7* is the MVT of X. 

k= t 

In particular, the essential MVT is exactly the union of the diagonal blocks in Figure 2. The 
parallelization of the case that we have just described can be easily done because it performs the 
same computations without the necessity of communication between the processors. In Figures 
1 and 3, the prior block may influence the posterior blocks for the MVT of X. For the Problem (jTJ 
of structurally nonsingular DAEs system, the fixed point iteration algorithm can give the same 
smallest dual-optimal pair with block fixed point iteration algorithm 113111 . 

Theorem 3.1. The above algorithm works correctly as specified and its complexity mainly de¬ 
pends on the block triangularization, finding the MVT and the block fixed point iteration as 
follows: 

(a) If X is irreducible, then the cost of the algorithm is 0(t + n(nlogn + 1 ) + qni 1 ) operations, 

(b) If X is reducible, then the cost of the algorithm is 0{p(NA> + N 2 logN ) + qN 2 ) operations, 
where q is the iteration number of fixed point for computing the global offsets c in X. 


Proof. The Correctness and termination of the algorithm follow from Lemma 3.6 in the com¬ 
panion paper [3JJ. 

Its complexity can be easily proved by Theorem 2.4 and Remark 2.5. □ 


Remark 3.2. In steps 3 and 4, for the existing MVT in X, we give the time complexity 0(pN 2 logN+ 
qN 2 ), which is superior to 0(pN 3 + qN 2 ) in M i In /fg[ I / ( AL the whole time complexity of their 
algorithms is 0(n 2 + nr + r(n + r)), where r is the number of adjusting the maximum weighted 
value from negative to zero to find a perfect matching in the bipartite graph. In particular, their 
index reduction techniques can be viewed as the behavior of signature matrix method by the bi¬ 
partite graph. However, they can not easily to achieve the block triangularization. Furthermore, 
Pryce only declared the Algorithm 3.1 infinitely many iterations in [21]. 


Example 3.1. Continue from Example 12. 71 we can mark the MVT with ’*’ and easily compute 
its gobal offsets by our extended 1,-method as follows: 



lil 

U\ 

d 

r 

r 

e 

z 

X 

C (D 

c < 2 ) 

c (3) 

c (4) 

c (5) 

/4 

rO* 



2 

0 




0 

0 

0 

0 

0 

h 


0 * 

2 


0 

0 



0 

0 

0 

0 

0 

fs 



0 * 

0 


0 


0 

2 

2 

2 

2 

2 

k 




0 * 


0 

0 


2 

2 

2 

2 

2 

h 





0 * 

0 


2 

0 

0 

0 

2 

2 

h 





0 

0 * 

2 


0 

0 

2 

2 

2 

h 







0 * 


2 

2 

0 

4 

4 

k 








0 * 

2 

2 

2 

0 

4 

d (D 

0 

0 

2 

2 

0 

0 

2 

2 






d< 2 > 

0 

0 

2 

2 

0 

0 

2 

2 






d (3) 

0 

0 

2 

2 

0 

2 

0 

2 






d (4) 

0 

0 

2 

2 

2 

2 

4 

0 






d (5) 

0 

0 

2 

2 

2 

2 

4 

4 
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We can get the offsets c = (2,2,0,0,2,2,4,4) and d = (4,4,2,2,2,2,0,0), for the variables in 
the original order given ©■ From step 5, the system Jacobian matrix J is 


m 2 

0 

0 

0 

r cos(9) 

sin(9 ) 

0 

° 

0 

m 2 

0 

0 

-Tsin(9) 

cos(9) 

0 

0 

0 

0 

Ml 

0 

0 

0 

-1 

0 

0 

0 

0 

J 

0 

0 

0 

C 3 

0 

0 

1 

sin{9) 

rcos{9) 

0 

0 

0 

0 

0 

0 

cos(9) 

-rsin(ff) 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 


Since det J = cos(9)(— tcos(9) 2 C 2 - Tsin(9) 2 Cf), whether this system is solvable or not depends 
on the angle 9, which is easy to check. Consequently, the extended 2,-method succeeds. 

4. Experimental results 

An efficient actual implementation of our algorithm is in Maple. The following examples run 
in the platform of Window and AMD Athlon(tm) IIX4-645 CPU 3.10GHz, 4.00GB RAM. First, 
we present an industrial application to illustrate our algorithm for the Five-axis Linkage CNC 
Machine YHV6025 models, see Figure 0 



There are 2446 differential and algebraic equations in the original physical system based on 
the Modelica language tool. We can get 245 differential and algebraic equations through the 
efficient elimination of trivial equations by means of symbolic simplification. The original sytem 
£ represents in Figure^ In Figure [7] we can obtain the block structure signature matrix for the 
corresponding £ based on our block triangularization algorithm. 
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Figure 6 : the original E 


Figure 7: the triangulated £ 


The triangulated signature matrix includes three blocks as follows. 
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^22 
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X 

X 

X 
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0 

X 
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y(3) 

^44 
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x 

x 

X 

X 

X 



where the blank in E, £ (1) , £ (2) and £ <3) means -oo, and the orders of signature matrices E* 1 ?, E^*, 
£®, E® and £^ 3) are 54, respectively. For each E (,) (i =1,2,3), we present the time for finding the 
global offsets (c) and (d) by the extended signature matrix method (ESMM), and then converting 
DAEs to the corresponding ODEs and computing the consistent initial values. In Table 1, we 
compare our method with signature matrix method (SMM) [2j] and weighted bipartite graph 
based on index reduction (WBGIR) |Q|. 


Table 1: Time for structural analysis of YHV6025 models 


Algorithm 

E * 11 



ESMM 

0.06036s 

0.05481s 

0.13926s 

SMM 

0.06942s 

0.06417s 

0.44486s 

WBGIR 

0.06823s 

0.06400s 

0.43745s 


From Table 1, our algorithm is more than three times faster than the other two algorithms for 
the large scale £ (3) . For the small and middle scale £ (1) and £ (2) , the advantage of our algorithm 
is not obvious compared with the other two algorithms. It shows that our algorithm is efficient 
for the large scale system of DAEs. 

We also give some results for large scale systems of differential algebraic equations via nu¬ 
merical simulation experiments. Assume that the original signature matrices E of systems have 
been preprocessed for the BTFs in Section lT2l and that the order of each diagonal block £,, in 
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the n x n matrix X is N in our experiments. In general, it is noted that the entries of X are the 
differential orders in DAEs. With loss of generality, we can randomly generate the sparse X for 
systems with BTFs properly as follows: 

• for each element in X 14 , randomly select an integer from {0, 1, 2, 3} according to its 
probability from {0.7, 0.1, 0.1, 0.1}, respectively; and X,,, = Xy,i = 2,3,...,/?; 

• for each element in Xyi, randomly select an integer from {- 1000 , 0 , 1 , 2 } according to 
its probability from {0.925, 0.025, 0.025, 0.025}, respectively; and X,,+i = X 12 , i = 
2,3,...,/?- 1; 

• the rest elements in X are -1000 (means - 00 ). 

We test the corresponding random trials for ESMM, SMM and WBGIR with N = {10,20,40} 
and n = 800 : 200 : 2400, and then calculate their constants in n • n v using the standard least- 
square method. The elapsed times of three algorithms are shown in Figure[ 8 ] respectively; some 
ratios of elapsed times are given in Figure [9] 





(c) WBGIR 


Figure 8: The elapsed times of structural analysis methods versus n. The constants in the legends 
are the fitted values /j in fin v for different N. 




(a) The ratio of elapsed times, i.e., 
versus n 


SMM’s elapsed time 
ESMM’s elapsed time 


,(b) The ratio of elapsed times, 
versus n 


i.e., 


WBGIR’s elapsed time 
ESMM’s elapsed time ’ 


Figure 9: Some ratios of elapsed times versus n. The constants in the legends are the fitted values 
fi in jun v for different N. 
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From Figure [8] we empirically know that the elapsed times of ESMM for the large scale 
nonlinear DAEs with different N are between O(n) and 0(n 15 ), and SMM seems to be 0(n 2 5 ). 
However, WBGIR seems to be 0(« 3 ). It is also noteworthy that WBGIR is very time consuming 
because of its recursive operations. Figure [9] shows that ESMM can reduce its elapsed time of 
sparse systems for fixed N by nearly 0(n ) (i.e., O(p)) compared to SMM and by 0(n 15 ) at least 
compared to WBGIR. In particular, the experimental results of ESMM are consistent with our 
complexity analysis in Theorem G. II 

Remark 4.1. For the large scale systems of sparse DAEs, we only consider the general reducible 
case for structural analysis. Moreover, in order to compare these three algorithms fairly, we 
construct the spare E of systems with BTFs properly. In practical applications, the signature 
matrix of sparse pattern S available from the large scale problems are often reducible. 


5. Conclusion 


In this paper, we propose an effective method to compute the large scale DAEs system in 
practice. We generalize the E-method with BTFs and combine with the block fixed point iteration 
for the general reducible DAEs systems. In particular, we exploit the shortest augmenting path 
algorithm for finding maximum value transversal, which can reduce the cost of computation. 
However, we apply the E-method as a basic tool, which relies heavily on square and sparsity 
structure. Thus it may be confronted with the drawback that can succeed yet not correctly in 
some DAEs arising from the specific types E l27ll . Some researchers also proposed some new 
computational techniques for high index nonlinear DAEs solving, such as the direct numerical 
computationlloj lhSl . index reduction [17, 25, 30h . and so on. 

In the future, we would like to consider the further research of hybrid symbolic and numerical 
computation with BTFs for large scale nonlinear systems of DAEs solving in practical industrial 
applicaitons, such as numerical algebraic geometry. 
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